"""
Created on 2018/2/3 18:09 星期六
@author: Matt  zhuhan1401@126.com
Description: 使用梯度上升法找到最佳参数
"""

from numpy import *
import matplotlib.pyplot as plt


def loadDataSet():
    dataMat = [];
    labelMat = []
    fr = open('testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2]))
    return dataMat, labelMat


def sigmoid(inX):
    return 1.0 / (1 + exp(-inX))


# 伪代码
# 重复R次：
#     计算整个数据集的梯度
#     使用alpha*gradient更新回归系数向量
#     返回回归系数
# dataMatIn：2维数组
def gradAscent(dataMatIn, classLabels):
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()  # 转置
    m, n = shape(dataMatrix)
    alpha = 0.001
    maxCycles = 500
    weights = ones((n, 1))
    for k in range(maxCycles):
        h = sigmoid((dataMatrix * weights))  # h为列向量
        # 计算真实类别和预测类别的差值,按照该差值的方向更新回归系数
        error = (labelMat - h)
        weights = weights + alpha * dataMatrix.transpose() * error
    return weights


# test One
# print("Test One")
# dataArr, labelMat = loadDataSet()
# weights=gradAscent(dataArr, labelMat)
# print(weights)


# 画出决策边界 即最佳拟合直线
def plotBestFit(wei):

    if(type(wei).__name__!='ndarray'):
        weights=wei.getA()
    else:
        weights=wei

    # weights = wei.getA()  # getA()用于将矩阵转换为数组  使用stocGradAscent0方法需要删除getA()
    dataMat, labelMat = loadDataSet()
    dataArr = array(dataMat)
    n = shape(dataArr)[0]
    xcord1 = [];ycord1 = []
    xcord2 = [];ycord2 = []
    for i in range(n):
        if int(labelMat[i]) == 1:
            xcord1.append(dataArr[i, 1])
            ycord1.append(dataArr[i, 2])
        else:
            xcord2.append(dataArr[i, 1])
            ycord2.append(dataArr[i, 2])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
    ax.scatter(xcord2, ycord2, s=30, c='green')
    x = arange(-3.0, 3.0, 0.1)
    y = (-weights[0] - weights[1] * x) / weights[2]
    ax.plot(x, y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

# test Two
# print("Test Two")
# dataArr, labelMat = loadDataSet()
# weights=gradAscent(dataArr, labelMat)
# plotBestFit(weights)


# 随机梯度上升算法
# 梯度上升算法每次更新回归系数时都要遍历整个数据集（“批处理”），如果样本过多，计算复杂度就太高了
# 改进方法：一次仅用一个样本点更新回归系数，称为随机梯度上升算法，这是一种“在线学习”算法
def stocGradAscent0(dataMatrix, classLabels):
    m, n = shape(dataMatrix)
    alpha = 0.01
    weights=ones(n)
    for i in range(m):
        # h 和 error都是相邻，上面的函数都是数值
        h = sigmoid(sum(dataMatrix[i]) * weights)
        error = classLabels[i] - h
        weights = weights + alpha * error * dataMatrix[i]
    return weights# 元组 tuple(3,)


# test Three
# print("Test Three")
# dataArr, labelMat = loadDataSet()
# weights=stocGradAscent0(dataArr, labelMat)
# plotBestFit(weights)

# 改进的随机梯度上升算法，使其能够更快地收敛
def stocGradAscent1(dataMatrix, classLabels, numIter=150):  # 默认迭代150次
    m, n = shape(dataMatrix)
    weights = ones(n)
    for j in range(numIter):
        dataIndex = list(range(m))  # bug 解决
        for i in range(m):
            alpha = 4 / (1.0 + j + i) + 0.01  # 每次迭代的时候调整α（吴恩达说不用。。。）,加入常数项使在多次迭代之后新数据仍有一定影响，α不会减小到0
            randIndex = int(random.uniform(0, len(dataIndex)))  # 随机选取更新
            h = sigmoid(sum(dataMatrix[randIndex] * weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            del (dataIndex[randIndex])
    return weights

# Test Four
# print("Test Four")
# dataArr, labelMat = loadDataSet()
# weights = stocGradAscent1(array(dataArr), labelMat)
# plotBestFit(weights)


# -------------下面进行病马死亡率的预测-----------

def classifyVector(inX, weights):
    prob = sigmoid(sum(inX * weights))
    if prob > 0.5:
        return 1.0
    else:
        return 0.0

def colicTest():
    frTrain = open('horseColicTraining.txt')
    frTest = open('horseColicTest.txt')
    trainingSet = [];
    trainingLabels = []
    for line in frTrain.readlines():
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(21):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[21]))
    trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 500)
    errorCount = 0;
    numTestVec = 0.0
    for line in frTest.readlines():
        numTestVec += 1.0
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(21):
            lineArr.append(float(currLine[i]))
        if int(classifyVector(array(lineArr), trainWeights)) != int(currLine[21]):
            errorCount += 1
    errorRate = (float(errorCount) / numTestVec)
    print('the error rate of this test is %f %%' % (errorRate * 100))
    return errorRate


def multiTest():
    numTests = 10;
    errorSum = 0.0
    for k in range(numTests):
        errorSum += colicTest()
    print('after %d iterations the average error rate is : %f %%' % (numTests, (errorSum* 100) / float(numTests) ))

# Test Five
# print("Test Five")
# colicTest()
multiTest()